% code used to create Fig. 4 (using CF=20000) Fig. S3
% in "An alternative pathway for delivery of 
% power by outer hair cells to cochlear traveling waves", Samaras and Meaud
%  plots the nonlinear response vs frequency for models M1, M2, M3, M4; 
% and plot experimental data from  Dewey et al. (2021)
% Julien Meaud
% julien.meaud@me.gatech.edu
% -------------------------------------------------------------------------

close all
clear all

global L

%% Colors
orange = [0.8500, 0.3250, 0.0980];
green = [140, 170, 0]./255;
yellow  = [0.9290, 0.6940, 0.1250];
purple = [0.4940, 0.1840, 0.5560];
grey = [0.6 0.6 0.6];
black = [0 0 0];
magenta = [1 0 1];
blue_CB = [26 133 255] ./ 255;
purple_CB = [212 17 89] ./ 255;
panelLabels = {'A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P'};
 
plotRL=0;
plotNormReCF=1;

% Patches for shading
xPatch = [0 55 55 0];
yPatch3 = [-10 -10 0 0];
yPatch4 = [0 0 2 2];
transparency_Patch = 0.5;

plotBMblack=0;
if plotBMblack 
    BMcolor = [0 0 0];
else
    BMcolor=blue_CB;
end
LS1 = '-';
LW=1.5;

leftX=0.12;
gapX=0.03;
height=0.25;
width=0.15;

nModel = 4;
r = 3;
c = nModel + 1;





Level=[30,50,70,80];

nLevel = length(Level);
transp = linspace(0.3,1,nLevel);

includeDCviscosity=2;

 deltaOHC = 0;
 CF=8000;  
 fSize=10;
 FreqLabel='Freq. re BF';

 for iModel=1:nModel
    for iLevel=1:length(Level)
    
        Level_curr=Level(iLevel);
   
        freqV=[2000:1000:30000]; 

        % load the data ---------------------------------------------------
        clear  uBM_M uDC_M uTMB_M
        for iFreq=1:length(freqV)
            freq=freqV(iFreq);
    
            if iModel==1  
                 load(sprintf('TimeDomain_CdcCohc/PTTD_G35_Act1_85B1_15A_rEps3_1_rDC_1000_Krl_BL_Cdc0_0_25_Cohc0_0_25_%gPT%g_SMALL',Level_curr,freq))           
            elseif iModel==2
                load(sprintf('TimeDomain_CdcCohc/PTTD_G35_Act3_15B1_15A_rEps3_1_rDC_5_Krl_BL_Cdc0_0_25_Cohc0_0_25_%gPT%g_SMALL',Level_curr,freq))                
            elseif iModel==3
                load(sprintf('TimeDomain_CdcCohc/PTTD_G35_Act1_65B1_15A_rEps3_1_rDC_1000_Krl_90BL_Cdc0_0_25_Cohc0_0_25_%gPT%g_SMALL',Level_curr,freq))                 
            elseif iModel==4
                load(sprintf('TimeDomain_CdcCohc/PTTD_G35_Act3_0B1_15A_rEps3_1_rDC_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_%gPT%g_SMALL',Level_curr,freq))                   
            end
        
            xMesh=ModelGeneratedVarsStruct.xMesh;
            freq=StimulusOptionsStruct.pureTone.F;
        
            uBM=ModelOutputsPostFFT.uBM;
            uDC=ModelOutputsPostFFT.uDC;  
            uTMB=ModelOutputsPostFFT.uTMB;
        
            uBM_M(:,:,iFreq)=uBM;
            uDC_M(:,:,iFreq)=uDC;
            uTMB_M(:,:,iFreq)=uTMB;
        
        end

        % find the best place using the low level data
        if iLevel==1
            indexFreq=find(freqV==CF);    
            [~,iBP(iModel)]=max(abs(uBM_M(:,2,indexFreq)));  
        end

        plotNormModel=CF;
        
        uBMreECP=squeeze(uBM_M(iBP(iModel),2,:)*1e7./convertSPLtoPressureInPascals(Level_curr)).*cos(deltaOHC);
        uDCreECP=squeeze(uDC_M(iBP(iModel),2,:)*1e7./convertSPLtoPressureInPascals(Level_curr));
        uTMBreECP=squeeze(uTMB_M(iBP(iModel),2,:)*1e7./convertSPLtoPressureInPascals(Level_curr));
        
        
        fig1=figure(1);
        set(fig1,'units','inches','position',[2 2 7 5.5])
        subplot('Position',[leftX+(width+gapX)*(iModel-1) 0.70 width height])
        plot(freqV./plotNormModel,20*log10(abs(uDCreECP)),'Color',[magenta transp(iLevel)],'LineWidth',LW,'Linestyle',LS1,'handlevisibility','off')
        hold on
        plot(freqV./plotNormModel,20*log10(abs(uBMreECP)),'Color',[BMcolor transp(iLevel)],'LineWidth',LW,'Linestyle',LS1,'Displayname',sprintf('%g dB',Level_curr))
          
        subplot('Position',[leftX+(width+gapX)*(iModel-1) 0.4 width height])
        plot(freqV./plotNormModel,20*log10(abs(squeeze(uBMreECP))),'Color',[BMcolor transp(iLevel)],'LineWidth',LW,'Linestyle',LS1,'Displayname',sprintf('%g dB',Level_curr))
        hold on
        plot(freqV./plotNormModel,20*log10(abs(squeeze(uTMBreECP))),'Color',[green transp(iLevel)],'LineWidth',LW,'Linestyle',LS1)
           
        subplot('Position',[leftX+(width+gapX)*(iModel-1) 0.12 width height])
        plot(freqV./plotNormModel,unwrap(angle(squeeze(uDC_M(iBP(iModel),2,:)./uBM_M(iBP(iModel),2,:))))/(2*pi),'Color',[magenta transp(iLevel)],'LineWidth',LW)
        hold on
        plot(freqV./plotNormModel,unwrap(angle(squeeze(uTMB_M(iBP(iModel),2,:)./uBM_M(iBP(iModel),2,:))))/(2*pi),'Color',[green transp(iLevel)],'LineWidth',LW,'Linestyle',LS1)   

    end
 end


% process the passive data (linear simulations) ---------------------------
for iModel=1:nModel

    if iModel == 1
        vars_Pass = load('Models_M1M2M3M4_FreqDomain/PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_1000_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat');
    elseif iModel == 2
        vars_Pass = load('Models_M1M2M3M4_FreqDomain/PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_5_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat');
    elseif iModel == 3
        vars_Pass = load('Models_M1M2M3M4_FreqDomain/PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_1000_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat');    
    elseif iModel == 4
        vars_Pass = load('Models_M1M2M3M4_FreqDomain/PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat');     
    end

    uBM_Pas = vars_Pass.ModelOutputsStruct.structural.uBM(iBP(iModel),:);
    uTMB_Pas = vars_Pass.ModelOutputsStruct.structural.uTMB(iBP(iModel),:);
    uTMR_Pas = vars_Pass.ModelOutputsStruct.structural.uTMR(iBP(iModel),:);
    freqPas = vars_Pass.StimulusOptionsStruct.freq;
   
    uBM_PasReECP = uBM_Pas.*1e7.*cos(deltaOHC)./convertSPLtoPressureInPascals(20);   
    uTM_PasReECP = uTMB_Pas.*1e7./convertSPLtoPressureInPascals(20);  
    uDC_Pas = vars_Pass.ModelOutputsStruct.structural.uDC(iBP(iModel),:);
    uDC_PasReECP = uDC_Pas.*1e7./convertSPLtoPressureInPascals(20);
    
    figure(1)
    subplot('Position',[leftX+(width+gapX)*(iModel-1) 0.70 width height])
    plot(freqPas./plotNormModel,20*log10(abs(uBM_PasReECP)),'Color',[BMcolor],'LineWidth',LW,'Linestyle',':','Displayname','Passive')
    hold on
     if iModel==1
        h2 = legend('AutoUpdate','off','Orientation','horizontal');
        h2.Location = 'northwest';
     end
     plot(freqPas./plotNormModel,20*log10(abs(uDC_PasReECP)),'Color',[magenta transp(iLevel)],'LineWidth',LW,'Linestyle',':')               
    
    if iModel==1      
         text(1.08,57,'BM','Color',BMcolor,'Fontsize',10)
         text(0.25,20,'OHC-DC','Color',magenta,'Fontsize',10)
    end
    if iModel==2
        text(7000./plotNormModel,20,'BM','Color',BMcolor,'Fontsize',10)
        text(5000./plotNormModel,60,'OHC-DC','Color',magenta,'Fontsize',10)
    end

    xlim([5 27.5]*1e3./20000)
    ylim([0 85])
    yticks([0:20:80])
    yticklabels([0:20:80])
    if ~plotNormReCF
        xticks([5:5:30])
        xticklabels([5:5:30])
    else
        xticks([0.25:0.25:1.75])
        xticklabels([]);
    end
    text(0.3,10,panelLabels{iModel},'Fontsize',16,'Fontweight','bold')
    yticks([0:20:80])
    if iModel == 1
        ylabel(['Disp. re ECP',char(10),'(dB re 1 nm/Pa)'])
        yticklabels([0:20:80])       
    else
        yticklabels([])
    end
    vline(CF./plotNormModel,'k--')
    set(gca,'Fontsize',fSize)

    if iModel==1
        title('M_1')
        h1=legend('boxoff');
        %h1.IconColumnWidth=10;
        h1.FontSize=10;

    elseif iModel==2
        title('M_2')
    elseif iModel==3
        title('M_3')     
    elseif iModel==4
        title('M_4')
    end

    subplot('Position',[leftX+(width+gapX)*(iModel-1) 0.40 width height])
    plot(freqPas./plotNormModel,20*log10(abs(uBM_PasReECP)),'Color',[BMcolor],'LineWidth',LW,'Linestyle',':','Displayname','Passive')
    hold on   
    if plotRL
        plot(freqPas./plotNormModel,20*log10(abs(uRL_PasReECP)),'Color',[purple_CB transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    else % plot TM
        plot(freqPas./plotNormModel,20*log10(abs(uTM_PasReECP)),'Color',[green transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    end

    xlim([5 27.5]*1e3./20000)
    ylim([0 85])
    yticks([0:20:80])
    yticklabels([0:20:80])
    if ~plotNormReCF
        xticks([5:5:30])
        xticklabels([5:5:30])
    else
        xticks([0.25:0.25:1.75])
        xticklabels([]);
    end
    text(0.3,10,panelLabels{iModel+c},'Fontsize',15,'Fontweight','bold')
    yticks([0:20:80])
    if iModel == 1
        ylabel(['Disp. re ECP',char(10),'(dB re 1 nm/Pa)'])
        yticklabels([0:20:80])
    else
        yticklabels([])
    end
    vline(CF./plotNormModel,'k--')
    set(gca,'Fontsize',fSize)
    if iModel==1
        text(7000./plotNormModel,50,'TM','Color',green,'Fontsize',10)
        text(15000./plotNormModel,22,'BM','Color',BMcolor,'Fontsize',10)
    end
    
    subplot('Position',[leftX+(width+gapX)*(iModel-1) 0.12 width height])
    plot(freqPas./plotNormModel,unwrap(angle(uDC_PasReECP./uBM_PasReECP))./(2*pi),'Color',[magenta transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    hold on
    plot(freqPas./plotNormModel,unwrap(angle(uTM_PasReECP./uBM_PasReECP))./(2*pi),'Color',[green  transp(iLevel)],'LineWidth',LW,'Linestyle',':')
    
    xlim([5 27.5]*1e3./20000)
    ylim([-0.5 0.5])      
    yticks([-1:0.25:1])
    if iModel==1
        yticklabels([-1:0.25:1])
    else           
      yticklabels([])
    end
        
    if ~plotNormReCF
        xticks([5:5:30])
        xticklabels([5:5:30])
    else
        xticks([0.25:0.25:1.75])
       xticklabels([0.25:0.25:1.75])            
    end

    text(0.3,-0.4,panelLabels{iModel+2*c},'Fontsize',16,'Fontweight','bold')
    vline(CF./plotNormModel,'k--')
    set(gca,'Fontsize',fSize)
    
    xlabel(FreqLabel)
    if iModel==1
        ylabel(['Phase re BM',char(10),'(cycles)'])
    end
        
           
end 

plotDewey_vs_freq  % plot experimental data+